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ABSTRACT 



A FORTRAN 77 computer code employing an adaptation of the 
finite differencing algorithm proposed by Brian was developed for the 
solution of transient heat conduction problems in cylindrical geome- 
tries. Validation of code was accomplished by comparison with an ana- 
lytic solution derived for a model with symmetric, linear boundary 
conditions. Accuracy of results for asymmetric and non-linear bound- 
ary conditions was determined by comparison with a similarly vali- 
dated code employing the explicit method. Code effectiveness was 
then demonstrated by conducting a transient temperature analysis for 
a simulated earth-orbiting satellite. Brian’s method demonstrated 
unconditional stability with associated significant reductions in execu- 
tion time compared to the explicit method. The effects of discretiza- 
tion error on the accuracy of results require further investigation. 
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I. INTRODUCTION 



A BACKGROUND 

Analysis of transient heat-conduction problems necessarily 
involves solution of the heat diffusion equation. In the absence of 
internal generation of energy, this is represented by 

v.(*vr)= P c,§- (11) 

Although analytic solutions may be obtained to certain simple transient 
problems [Ref. l:p. 212], many of the practical problems encountered 
by engineers involve one or a combination of nonlinearities, complex 
geometries, complex boundary conditions, or systems of coupled par- 
tial differential equations of which any one may necessitate the use of 
numerical methods. Although Monte-Carlo (probability sampling) and 
finite-element methods have been applied to the solution of heat-con- 
duction problems, the approximation of partial derivatives by finite 
differences presents a straightforward and popular approach [Ref. 2:p. 
471]. Two fundamental finite-difference techniques are the implicit 
and explicit methods and their derivatives. 

B. PROBLEM 

The analysis of transient heat-conduction problems in three 
dimensions is of great importance to many fields of science and engi- 
neering, but solutions are often quite costly in terms of computational 
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effort and time consumed in their solution. Numerical solutions to 
these problems frequently tax the computational and storage capaci- 
ties of available computers [Ref. 3:p. 367]. 

The explicit method offers the benefits of relatively low storage 
requirements and simplicity in formulation, coding, and execution. 
However, it is only conditionally stable and is subject to a limitation on 
the maximum timestep which may be employed. If the maximum 
timestep is exceeded, the solution experiences numerically induced 
oscillations which may cause the solution to diverge from the correct 
result [Ref. l:p. 214]. This timestep restriction causes the explicit 
method to be inefficient in terms of computational effort for the analy- 
sis of lengthy transient periods. 

The implicit method offers the advantage of being unconditionally 
stable and therefore exempt from timestep restrictions. However, this 
method calculates node temperatures by the simultaneous solution of 
the nodal heat balance equations. Thus, savings gained from the lack of 
a timestep restriction are often quickly offset by the increased 
computational effort associated with the required matrix inversions. 

The implicit alternating direction method (ADI) is an uncondi- 
tionally stable extension of the implicit method which reduces the 
number of required computations by transforming the problem into 
tridiagonal form. However, a direct extension of this method into 
three dimensions becomes only conditionally stable, is restricted to a 
maximum timestep in a fashion similar to the explicit method, and 
demonstrates similar shortcomings [Ref. 4:p. 453]. 
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C REQUIREMENT 



A finite difference technique which is unconditionally stable in 
three dimensions and unencumbered by the requirement to manipu- 
late large matrices is required for the solution of transient heat- 
conduction problems in three dimensions. The algorithm proposed by 
Brian [Ref. 3] offers such a technique. 



D. OBJECTIVES 

1. Adapt the algorithm proposed by Brian to the solution of three- 
dimensional transient heat-conduction problems in cylindrical 
geometries. 

2. Evaluate the relative accuracy of the results obtained using Brian’s 
Method by comparing them against those obtained using the 
explicit method. 

3. Evaluate the speed of execution of Brian’s Method relative to that 
of the explicit method in terms of central processing unit (CPU) 
seconds. 

4. Illustrate the relative benefits derived from the Brian’s Method by 
applying it to the solution of an application problem involving 
long-period transient behavior. 
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H. THEORY 



A FUNDAMENTALS 

1. Material Properties 

The capacity for energy storage and rate of energy transfer 
within a medium are functions of the material’s thermophysical prop- 
erties. Isotropic materials, such as metals and alloys, are uniform in 
nature and their properties generally demonstrate only negligible 
directional and temperature dependence. However, layered or strati- 
fied materials, such as glass-reinforced plastic (GRP) or laminates, may 
demonstrate strong dependence on both direction and temperature. 
While directional dependence of properties is generally linear in 
nature, temperature dependence is strongly nonlinear and will not be 
addressed. 

a Heat Capacity 

The rate of change of energy storage per unit volume 
within a solid, isotropic medium is determined by the product of its 
density (Kg/m 3 ), specific heat (KJ/KgK), and the time-rate of change 
of its temperature : 

E '- = P C "W (21) 

Both density and specific heat may be temperature and/or direction- 
ally dependent properties. 
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b. Thermal Conductivity 

The thermal conductivity (W/mK) is a transport property 
(i.e., one affecting the movement of energy) of the material and is 
defined in the x-coordinate direction as 







( 2 . 2 ) 



Isotropic materials generally demonstrate relatively constant thermal 
conductivity over broad temperature ranges, whereas metallic lami- 
nates demonstrate a strong directional dependence within the plane 
of each different material. Composite materials show both strong 
directional and temperature dependence. Temperature-dependent 
material properties will not be addressed. 

c. Thermal Diffusivity 

The thermal diffusivity (m 2 /s) is a grouping of the mate- 
rial density, specific heat, and thermal conductivity and represents 
the controlling transport property for transient heat conduction [Ref. 
l:p. 41] 



a = 



k 

pc P 



(2.3) 



d. Convection Heat Transfer Coefficient 

The convection heat transfer coefficient (W/m 2 K) 
encompasses all of the effects that influence the convection mode and 
is dependent upon boundary conditions, surface geometry, the nature 
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of the fluid motion, and a number of other fluid thermodynamic and 
transport properties (Ref. l:p. 8]. 

e. Properties Affecting Radiation 

Various surface properties affect the radiative character- 
istics of a solid material. Directionally integrated or hemispherical val- 
ues are often adequate although these properties are still strongly 
dependent upon the wavelength of the radiant energy. Where appro- 
priate, consideration of the wavelength dependence will be indicated 
by a subscripted Greek letter lamda (X). 

(1) Emissivity. The emissivity (e) is a surface property 
representing the ratio of total energy emitted by an actual surface and 
an identical black surface at the same temperature. The emissivity may 
be altered by one or a combination of mechanical processes (e.g., 
sandblasting or metal-spraying) or by application of various films and/ 
or coatings, either by design or through service conditions. 

(2) Surface Properties. Absorptivity (a) and reflectivity 
(p) are surface properties representing the fractions of total incident 
energy which are respectively absorbed and reflected by the surface. 
For opaque surfaces, the sum of the absorptivity and the reflectivity 
equals one 



a+ p= 1 (2.4) 

In addition to absorption and reflection, translucent materials allow a 
fraction of incident radiation to be transmitted through the material. 
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Consequently, an additional property, transmissivity ft), must be 
incorporated into Equation 2.4, yielding 



a+ p+ t = 1 (2.5) 

(3) Radiation Heat Transfer Coefficient. The radiation 
heat transfer coefficient (W/m 2 K 4 ) is a collection of terms resulting 
from the linearization of the radiation heat transfer calculation from a 
gray, diffuse surface in extensive surroundings and is represented as 



h r = e a (T s + T sur )(T 2 s + rL) (2.6) 

where o is the Stefan-Boltzman Constant (5.67 x 10' 8 W/m 2 K) and T s 
and T sur are the absolute temperatures of the surface and the sur- 
roundings, respectively. As is apparent from Equation 2.6, this term is 
strongly temperature dependent and non-linear in nature. 

2. Heat Transfer Rate 

The heat flux (W/m 2 ) is rate of heat transfer per unit surface 
area and may be represented by 



q 



// 



E_ 

A 



(2.7) 



The heat transfer rate (W/m 2 s) is thus the product of the heat flux and 
the surface area normal to the direction of heat flow [Ref. l:p. 4]. 
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3. Rate Equations 

The three principle modes of heat transfer are conduction, 
convection, and radiation. The rate of heat transfer attributed to each 
mode is described by the appropriate rate equation. 

a Conduction Equation 

The governing equation for conduction is Fourier’s Law, 
given for the x-direction by 



Q 



// 

X 




(2.8) 



b. Convection Equation 

The relation governing convection is known as Newton’s 
Law of Cooling and is given by by 



q"c= h c (T s -T sur ) (2.9) 

where T s is the surface temperature and T^ is the temperature of the 
surrounding fluid. 

c. Radiation Equation 

The maximum flux at which radiation may be emitted by 
a surface to an infinite, black surrounding is determined by the Stefan- 
Boltzman Law 



q" r = K (T s - rj (2.10) 

where T s is the surface temperature and Tsur is the temperature of 
the surrounding enclosure [Ref. l:p. 9]. In reality, objects exchange 
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energy with surrounding enclosure boundaries in accordance with the 
relation 

q"=^F.A.h r (T s - T sur ) (2.11) 

.■ ' • 

where F is the view factor representing the fraction of the radiation 
emitted by the surface of interest that is intercepted by the i th surface 
of the confining enclosure. Implicit in Equation 2.11 is the assumption 
that all surfaces of the surrounding enclosures are black. 

4. Heat Balance — The First Law of Thermodynamics 

The temperature at a given point may be determined from 
establishment of an arbitrary control volume around the point (Fig- 
ure 1) followed by application of the first law of thermodynamics (the 
law of conservation of energy). Excluding heat which is generated 
within the confines of the control volume (internal generation), this 
law may be simply stated as 

The rate at which thermal and mechanical energy enters a control 
volume minus the rate at which this energy leaves the control vol- 
ume must equal the rate at which this energy is stored in the 
control volume. [Ref. l:p. 12] 

With the addition of any heat genetrated internal to the control vol- 
ume, the first law may represented in the form of the heat balance 
equation: 



E = E. 

& in 



E - 

^ gen 



( 2 . 12 ) 
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Figure 1 . Control Volume 



where 

E a = Rate of energy storage 
« = Rate of energy gain 

P 

° ui = Rate of energy loss 

r 

*" = Rate of energy generation. 

R FINITE DIFFERENCE APPROXIMATIONS 
1. General Approach 

The numerical solution of transient, heat-conduction prob- 
lems by digital computer requires the transformation of the heat con- 
duction equation into a form which allows differentiation to be 
approximated by finite differences [Ref, 2:p. 472]. This may be accom- 
plished by dividing the problem into a number of discrete points, or 
nodes, each surrounded by an associated control volume (Figure 2(A)). 
Conducting an energy balance on each control volume yields a system 
of linear algebraic equations which describes the temperature field. In 
the absence of round-off error, the solutions obtained from finite 
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difference approximations tend to the analytic solution as the step size 
and grid spacings tend to zero [Ref. 4:p. 449]. 




Type 1 — End-layer, inner-surface 
Type 2 — End-layer, mid-radius 
Type 3 — End-layer, outer-surface 
Type 4 — Mid-layer, inner-surface 
Type 5 — Mid-layer, mid-radius 
Type 6 — Mid-layer, outer-surface 

Figure 2. (A) Reference Grid; (B) Node Types 
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cl Geometry 

Figure 2(B) depicts a right-circular cylinder to which has 
been applied a polar-cylindrical coordinate system and an arbitrary 
grid which divides the cylinder into a number of nodes, each sur- 
rounded by an associated control volume. Nodes may be characterized 
as belonging to one of six types and are located at the grid intersection 
points. Each node may be associated with a specified layer, ring, and 
ray. 

(1) Layers. Layers lie parallel to the x-y plane and are 
numbered sequentially beginning in the x-y plane and proceeding in 
the direction of the positive z-axis. Layers are spaced a distance Az 
apart. Note that the end layers are in the plane containing the ends of 
the cylinders. 

(2) Rings. Each layer is divided into an arbitrary number 
of concentric rings spaced a distance Ar apart. The node at the three 
o’clock position when viewed from above has been arbitrarily desig- 
nated as node one. Other nodes are numbered sequentially in a 
counter-clockwise direction from node one. 

(3) Rays. Rays emanate outward from the origin in the 
plane of the layers. Nodes lie along the rays at intervals of Ar and are 
numbered sequentially from the inner surface outward. 

b. Subscripting Convention 

To simplify notation in the development of equations, the 
designation of nodal grid coordinates will follow the convention con- 
tained in Table 1. 



12 



TABLE 1 



NODAL SUBSCRIPT CONVENTION 



Subscript 


Location 


Grid Coordinates 


N 


Northern (upper) neighbor 


[r,<Mz+l)] 


S 


Southern (lower) neighbor 


[r,<(),(z-l)] 


E 


Eastern (right) neighbor 


[r,(<)>-l),z] 


W 


Western (left) neighbor 


[r,(4>+l),z] 


I 


Inner neighbor 


[(r- l),4>,z] 


O 


Outer neighbor 


[(r+l),<{>,z] 


e 


Edge 


(varies) 



* All grid coordinates are relative to the node of interest when viewed 
from the inside of the cylinder. 



In view of the inverse symmetric relationship existing 
between equivalent nodes at the upper and lower surfaces, the sub- 
script e will be used to designate temperatures, fluxes, and properties 
associated with the end layers, thus reducing the number of required 
equations by half. The reader must substitute the appropriate subscript 
(i.e., N or S) depending upon the location or application of the equa- 
tion. This issue will be further elaborated upon in the equation devel- 
opment example. 

c. Superscripting Convention 

Superscripts, as contained in Table 2, will be employed 
to indicate at which timestep a value is considered. 
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TABLE 2 



NODAL SUPERSCRIPT CONVENTION 



Superscript 


Timestep 


Time 


n 


Current 


t 


n+1 


Next 


t + At 


n+1/2 


Intermediate 


t + At/2 


<L Example 







Consider the one-dimensional heat conduction problem 
represented in Figure 3 wherein T E = T w , heat is being conducted in 

the radial direction only, and the control volume is of unit thickness. 



T 

A o 




Figure 3. One -Dimensional Heat Conduction Problem 

The appropriate heat balance equation is 

E a =E,-E 0 (2.13) 
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The time-dependent storage term ( E a ) may be represented by its 
finite difference approximation 




(2.14) 



where At is the step size, T 11 is the node temperature at the old (n) 
timestep, T n+1 is the node temperature at the new (n+1) timestep, 
and the volume is approximated as 



The heat rate in the radial direcrion at the inner bound- 
ary ( E I ) of the control volume may be approximated as 



is the area of the control volume in the direction of heat flow. Simi- 
larly, the heat rate at the outer boundary may be approximated as 



E,= - k r A 



(T* -T n+1 ) 



(2.16) 



where 



A t = (r- -^f) A0 



(2.17) 




(r;- r" +1 ) 



(2.18) 



A r 



where 



A 0 =(r+-^)ArA<p 



(2.19) 
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Substituting these into Equation 2.13 yields the finite 
difference form of the heat equation 



A r. 



pC p rA rA(p (T" +1 - T") M r “ — ^ 



Ar A r 

k r (r+ A r)A(p 



- r n ) 



( 2 . 20 ) 



A r 



(r;- r ) 



The steps leading up to and including the formation of Equation 2.20 
are common to both the explicit and Brian’s methods. 

Similar equations may be obtained for other nodes and 
boundary conditions. 

2. The Explicit Method 

The explicit method is a direct application of the general 
approach to finite difference analysis detailed in the preceding sec- 
tion. Equations describing the temperature at the new time level are 
formed for each node by conducting a heat balance on the associated 
control volume. Temperatures at the new time level may then be 
approximated explicitly by substituting temperatures at the old time 
level into the appropriate node equation. The explicit form of the one- 
dimensional heat conduction problem of the preceding section may be 
obtained from Equation 2.20 by dividing by the lead coefficient on the 
left-hand side, substituting for the grid Fourier number, and solving 
for the new temperature, yielding 



rji +1 r-p, I 



+ Fo 



' , Ar, 

(/--—) 

A r 



(T, n - T n )+ 



t . A r , 
(r + — ) 

A r 



( Tq T ") 



( 2 . 21 ) 
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where 




( 2 . 22 ) 



Defining the constants 




(2.23) 



and 



c. 



2 



A r 



(2.24) 



and grouping the T 11 terms outside of the the parentheses yields 



from which the stability criterion on the timestep may be determined 
from the relation 



[Ref. l:p. 215] 

3. Brian’s Method 

In Brian’s Method, the node types are the same as depicted 
in Figure 2. However, the temperature at the new time level is formed 
from a combination of values taken from the old time level and three 



r +1 = Fo(cj; + C 2 T ;) + [l - Fo( Cl + C 2 )]T 



(2.25) 



[1 - Fo( Cl + c 2 )] > 0 



(2.26) 
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directionally dependent temperature arrays formed at an intermediate 



time 



t' = t + (2.27) 

Although each of these arrays requires a full complement of supporting 
node equations, the equations are quite similar and enjoy a common 
derivation with the explicit method for the initial stages of their 
development. 

Each intermediate array is associated with one arbitrarily 
assigned coordinate direction and has its temperatures denoted by the 
superscripting of one, two, or three asterisks. Calculation precedence 
must follow the order of star level, two-star level, three-star level, and 
new temperature arrays! 

Within each intermediate calculation process, temperatures 
are determined by solution of a tridiagonal matrix of coefficients of 
nodes lying along the specified coordinate axis. For example, equations 
for the r-direction nodes at the star level in a one-dimensional heat 
conduction problem are of the form 

ATj+ BT’+ CT' 0 = D(T n ) (2.23) 

Although Brian’s method is conceptually no more difficult 
than the explicit method, it involves considerably more algebra and is 
best illustrated in the form of the code development model contained 
in the following chapter. 
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III. DERIVATION OF EQUATIONS 



A MODEL FOR CODE DEVELOPMENT 

1. General 

In this chapter, detailed descriptions of the explicit and 
Brian’s methods will be presented in the form of the model employed 
for the development of code. 

2. Problem Description 

Selection of a model with simplicity and symmetry of geome- 
try and boundary conditions was made to facilitate error checking and 
the comparison of results during code development. The preliminary 
model consisted of a right-circular cylinder of uniform properties with 
adiabatic ends and specified initial temperature distribution (initial 
condition) and constant inner-surface temperature of zero degrees 
Celsius. An arbitrary grid consisting of five layers in the z-direction, 
each with three and four nodes respectively in the radial and phi- 
directions, was applied to the cylinder in a fashion similar to that 
depicted in Figure 4. At time t = 0, a uniform, constant heat flux was 
applied to the outer surface. Subsequent modifications to the model 
added radiative and convective boundary conditions at all surfaces as 
well as directionally-dependent thermal conductivities. Temperatures 
are in terms of the rise above the initial temperature (excess 
temperature) [Ref. l:p. 100]. 
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Specified Temperature 



Adiabatic End 




Adiabatic End 



Figure 4. Model for Code Development 

R EXPLICIT METHOD EQUATIONS 

The process for development of equations using the explicit 
method is similar for all node types. For illustrative purposes, we will 
consider a Type 1 node located at the inner edge of the bottom sur- 
face of the cylinder as depicted in Figure 5. Complete equations for all 
node types are included as Appendix A to this report. Directionally 
dependent thermal conductivities and radiative and convective bound- 
ary conditions are included in this development. Heat flow into and 
out of the associated control volume is represented in Figure 6. 
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(r + A r),(p, z 




r,{(p + A 0), z r,(<p - A <p),z 



A. Plan View 



r,0 , ( z + Az ) 



r,(0 + A0 ), z 



-• — 
r,<p, z 



— • 

r,(0 - A0), z 



B. Frontal View 

Figure 5. Bottom Layer, Inner Surface Node Configuration 
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E 






E, 




E, 



o 



Figure 6. Heat Balance 

The heat balance equation is 



E * = E , ~ E n ~ E s ~ E e ~ E w ~ E o 



( 3 . 1 ) 



Application of appropriate finite difference approximations results in 

+ iT n ) 



n + l _ n (7* - r") w 

pcy- a » ( 3 . 2 ) 



At 



Az 



* * 



r A(p 



(T 0 n - T n ) 

+ krAr ~ 9 ~A~r E s + E ‘ 



where 



V = 



( r+ ^~)A rA(p A z 



( 3 . 3 ) 
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( r + 4p)A0 Az 



o 



2 



(3.4) 



A rAz 



(3.5) 



A 



(r+^)ArA<p 



2 



(3.6) 



We note in Equation 3.2 that the conduction terms are evaluated 
at time level n in order to allow explicit solution of the temperatures. 
Evaluation of E e and E l reveals that these terms are comprised of the 
sum of the incident surface heat flux and the terms associated with 
heat loss (or gain) by convection and radiation. The convection term is 
linear and may be approximated as 



where h c is the convection coefficient at the appropriate surface and 
Too is the temperature of the surrounding fluid. The radiation term is 
non-linear and cannot be directly applied as a function of T n+1 . A lin- 
earization of the radiation equation may be accomplished by defining 
the radiation coefficient at the previous (n) timestep 




(3.7) 



h r = ecj(T n + T sur )[{T n f + TL)] 



(3.8) 



where 




(3.9) 
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which reduces the approximation of the radiation term to the first 
power of T n+1 and T sur is the temperature of the surrounding envi- 
ronment to which it radiates. Thus, the boundary conditions may be 
represented as 



E s = A,[q" s - h u (T"‘ - T.) - h, s (T”' - r„)] 



(3.10) 



and 



E] = A,[q",-h'(T-"- TJ- h r (T'"-T m )] ( 3 . n) 



where 



rAtpA z 



(3.12) 



Substituting for E s and E t into Equation 3.2 yields 



pc P v 



^ rj~> ^ +1 f 

A 1 



= k,A 



(J N - T ) 
A z 



4 - k, A. 



(J 



+ T 

E W 



-2 T n ) 



rA<p 



(3.13) 



(T" - r") 

+ k,A 0 -\- r + A^q"; + h c {T„- T n+l )+ h rs {T sur - T" +1 )] 

+ A,[q" l+ h c {T„- T n+i )+ h r (T sur - T" +1 )J 
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Substituting for the volume and area terms yields 



pc, 



(r + -4p-)ArA0 Az ( 7" +1 _ 7 ") k 2 (r + ^-)ArA<p ( 7 ^ _ 7 ") 



Ar 



Az 



(3.14) 



k^ArAz (T n E +T n w -2T n ) k,(r + Af)A4> Az ( 7 ^- 7 ") 



rA0 



Ar 



. Ar . . . , 

(r + -r-)ArA(l> 

+ W' s + h c (T„-T n+1 ) + h r (T sur - 7 n+1 )] 

2 s s 



+ rA 2 Az [ + K ' <r_ -r ■*') + a, ; (r,„ - r *')] 



By defining the directional thermal diffusivity 



a 2 = 



pc, 



(3.15) 



and grid Fourier number 



Fo = 



a Ar 
A z 2 



(3.16) 



we may simplify by arbitrarily dividing by the lead coefficient of the z- 
direction difference approximation to obtain 
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k^Az 2 



1 

2 Fo 






- T ") = (T n N - T ") + 



2k z r (r + ^-)M> ~ 



(T 



l+K-in 



k r (r + -y-)A z 2 

+ b — v n o-T n ) 

k t (r + ^p)A r 2 



+ -jZ-W' s +h c ( T„-T n+l ) + h r (r w -f +1 )] 

k 2 4 s s 



rAz 2 



k 2 (r + —)Ar 



~rz [<f, + /i e (T„-T n+ 1 )+h r (T sur -T n+ ')] 

A r ^ . _ / / 



Further simplification may be realized by defining the constants 

£ A z 2 

2 ''( r+ ^f-) A( A 2 



£,( r + "4p)Az 2 

^,(r + -^p)A z 2 



C 3 = 



Az 



C 4 = 



rA z 2 



, . Ar. . 
^ 2 (r+ — )A r 



(3.17) 



(3.18) 

(3.19) 

(3.19) 

(3.20) 
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Substituting into Equation 3.16 and rearranging to yield 



2 Fo 



(t - t ") = ( r ; - r n ) + C] (t; +T ;-2T n ) + c 2 (t ; - r > 



(3.21) 



+ c 3 [q " s - h c (T„-T n+l )-h r (T^-T^)} 

3 3 s s 

+ c 4 [ q"i — T n+1 )] 

Expansion and regrouping of terms yields the final form of the 
node equation for the temperature at the (n+1) timestep 

r n+1 = ~}r{ T " + 2 Fo[(T* - T n )+ c x (T n E + r; -20 + c 2 (r o n - r") (3.22) 

+ c 3 ( q" s + h c t x + h r rj 

3 s s 



+ c 4 (Q'i + h c T„+ h r T j- 



where 



c 5 = l + 2Fo[c 3 (h c s + /i r H c 4 (h c ' + h r ) ] (3.23) 

Thus, the temperature for each node at the successive time step 
may be solved for explicitly in terms of node temperatures at the pre- 
vious timestep. 

C BRIAN’S METHOD 

1. Step One (z-Direction) 

This section addresses procedures governing the develop- 
ment of equations associated with the first intermediate, or star-level , 
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temperature array to which the z-coordinate direction has been arbi- 
trarily assigned. Equation 3.21 may be applied to step one by substi- 
tuting asterisk-superscripted variables for all n+1 superscripted 
variables and those n-superscripted variables which lie along the z- 
coordinate direction. Because the intermediate time step is half the 
overall value, the effective grid Fourier Number is half that employed 
in the explicit method. Thus, for step one of Brian’s Method, Equation 
3.21 becomes 



Expanding and grouping like terms of starred quantities on the left- 
hand side of the equation yields 




+ c,[q" s - h e {T-TJ-h r {T'- rj ] 



+ c A [q",- h c ( r-TJ-h r (7*- T S J] 



AT) + BT + cr N = D 



(3.25) 



where 



A= 0 



(3.26) 




(3.27) 



C = - Fo 



(3.28) 
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D = T" + Fo[ Cl (T n E + T n w - 27’")+ c 2 (T„ - T*) 



(3.29) 



+ c ( q" s + h c T„+ h r T sur ) 

3 s s 

+ c A {q", + h c T„+ h r T W )J 

This process is repeated until similar equations are derived for each 
node type. 

2. Step Two (^-Direction) 

Derivation of step two equations is an extension of the proce- 
dures employed in the preceding section with the phi-axis now arbi- 
trarily designated the two-star coordinate direction. In this step, 
beginning with Equation 3.24, which already contains the initial step 
one substitutions, we add further modifications by substituting two- 
star subscripted variables in the same manner as before. This yields 

i(T"-T ') = (?■;- n + q(r"+ + c jt;~ t") 0.30) 

+ c,[V s - K JJ" - T.) - K (T - - r„)] 

+ h c (r-- T.)-h,(T"- r,„)] 

Again, expanding and grouping like terms of two-starred quantities on 
the left-hand side of the equation yields 

AT” + BT” + CT” = D (3.31) 
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where 



A=-c x to (3. 31 ) 

5=1+ Fo^2c l + c 2 ( h c + h r ^) + c 3 ( h c , + h r ) J (3.32) 

C = -c,Fo (3.33) 

d = r n + Fo[(t* n - r ) + c 2 (r 0 " - r") (3.34) 

+ c [ q '' + h c T„ + h r T sur ] 

* * S 5 

+ c 4 (<7 ” + h e T„+ h r T W )J 

This process is again repeated until similar equations are derived for 
each node type. 

3. Step Three (r-Direction) 

Derivation of step three equations is an identical extension of 
the procedures employed in the preceding section, with the r-axis 
now arbitrarily designated the three-star coordinate direction. In this 
step, we begin with Equation 3.30 and further modify by substituting 
three-star subscripted variables. This yields 
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r") = (r; - o + q(r” + t“-2T”) + c 2 (7”' - t***) (3.35) 
+ c 3 [<?" - h c {T~-Tj-h r (r**- rj] 

e < 

+ c 4 [<?" - h e {T'“- TJ- h r (T'“- rj] 

Again, expanding and grouping like terms of two-starred quantities on 
the left-hand side of the equation yields 

at;** + 57**’ + cr;** = d (3.36) 

where 

A= 0 (3.37) 

5=1+ Fo[c 2 + c 2 {h c ^ + h r )+ c 4 ( h c ^ + h r (3.38) 

C = ~c 2 Fo (3.39) 



D= T n + Fo[(T‘ n - r‘)+ c,(r;*+ T” -2T”) (3.40) 

+ c, ( + h c T„+ h r T sur ) 
s s 

+ c 4 ( < ? / / + A. /-+ ^r / 7’ Jur )j 

This process is repeated until similar equations are derived for each 
node type. 
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4. Step Four (New Temperatures) 

Step four is a direct extension of the procedures employed in 
the explicit method. In this step, the n+1 timestep equations are 
derived from the explicit equations by substituting starred values for 
the associated node temperatures for node under consideration. Con- 
tinuing with our example, we substitute into Equation 3.22 to obtain 

T H+l = T~ 5 {T n + 2 Fo[(T‘ n - T‘) + Cj(r‘* + r;*-2T")+ c 2 (T' 0 “- T*”) (3.41) 

+ c (q" + h c L+ h r T sur ) 

J 15 s s 

+ c 4 (<7 ", - h c T„- hrT^yj 

where 

c 5 = 1 + 2 Fo[c ( h c + h r )+c(h c + h r )] (3.42) 

s s i i 

Step four equations for each of the six types of nodes are 
derived in an identical fashion. 

Complete node equations in the polar-cylindrical form of 
Brian’s method are included in Appendix B for each of the six node 
types identified in Figure 2. 



32 



IV. CODE DEVELOPMENT 



A ANALYTIC SOLUTION 

Analytic solutions to the model described in Chapter III, para- 
graph A.2, were obtained for the range of zero to 300 seconds in five- 
second increments. The calculations of eigenvalues and the resulting 
solutions were made using MATHCAD V2.0 on an INTEL 80386-based 
personal computer. The calculated eigenvalues were verified with tab- 
ular results [Ref. 5] and those obtained by an appropriate IMSL 
subroutine [Refs. 6-10]. Equations describing the complete solution 
are contained in Appendix C. 

R EXPLICIT METHOD 

1. Requirements 

The coding of the explicit method is relatively straightfor- 
ward and subject to only two notable requirements: 

• That the selected timestep must not exceed the theoretical 
maximum, and 

• That each node be considered for temperature calculation once 
per timestep. 

2. General Procedure 

a Timestep Restriction 

The maximum allowable timestep varies with control 
volume dimensions. Exceeding the allowable timestep for any one 
node will result in the eventual loss of stability of the entire system. 
The rapidity with which this occurs is a function of the affected node’s 
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significance in the path of the dominating heat flux. Consequently, a 
maximum timestep must be determined for each node type (Figure 2) 
and the most restrictive taken as the limiting value. This was effected 
through use of the subroutine STABIL, which calculated maximum 
allowable timestep for each node-type and returned the most restric- 
tive as the argument DTMAX. A listing of equations defining stability 
requirements for each node type is contained as Appendix D. 
b. Indexing 

The method by which incrementation through the grid 
system is implemented is entirely arbitrary as long as the temperature 
is calculated once per timestep for each node on the grid. In systems 
with a dominant heat flux, convergence will occur more quickly if the 
dominant direction is incremented first. Appendix E contains the 
explicit method code as developed for this study. A polar-cylindrical 
coordinate system was arbitrarily applied with the z-axis aligned with 
the central axis of the cylinder. Indexing commenced at the inner 
edge of the bottom surface and proceeded up the z-coordinate to the 
upper surface. The grid was then indexed in <f>, followed by r, until the 
entire grid had been covered, at which time the old temperature 
matrix was replaced by the new temperature matrix, the time counter 
indexed by At, and the process repeated until the end of the desired 
period of evaluation was reached. 
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C BRIAN S METHOD 
1. Requirements 
a Times tep 

The timestep may be selected as befits the duration of 
the analysis and the expected temperature gradients. However, the 
discretization error is a function of the square of the timestep. On the 
other hand, the property of unconditional stability allows the flexibility 
to initially choose a large timestep for macro-analysis followed by 
application of a smaller timestep for a detailed look at regions of 
interest (e.g., an area where a predefined threshold is exceeded). 
b. Indexing 

A separate temperature matrix must be developed for 
each coordinate direction at the level of the intermediate timestep. 
Due to this directional nature, the indexing requirements are more 
complex and restrictive. Once assignment of a coordinate axis to an 
intermediate matrix (e.g., star-level) has been made, the convention 
must be adhered to for the duration of the timestep. One acceptable 
method for indexing is outlined below and the companion code is 
included as the STEP(#) subroutines to the Brian’s Method code con- 
tained in Appendix F. 

(1) Star Level (z-Cirection). The z-coordinate direction 
was arbitrarily assigned to the star-level intermediate matrix. 
Consequently, equations were developed for nodes lying along the z- 
axis in the fashion described in Section C of Chapter III. Appropriate 
indexing was accomplished by designating an arbitrary node on the 
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inner ring of the bottom surface of the cylinder (Figure 7) as a node of 
reference (node 1). An equation similar to Equation 3.25 was devel- 
oped for this first node and the values of the A, B, C, and D coefficients 
calculated and stored in four associated vectors. The pointer was then 
indexed in the z-direction to select the node immediately above 
node 1 (i.e., the northern neighbor when developing the equation for 
node 1). The process was repeated for the second node and all other 
nodes lying on the selected z-axis 




Figure 7. Star-Level Indexing Direction 



When the values of the coefficients for the last node had been calcu- 
lated and stored in the appropriate vectors, a tridiagonal matrix solver 
was called to simultaneously solve the system of equations. The resul- 
tant node temperatures were then saved in the star-level intermediate 
matrix and the. pointer indexed to the next adjacent node in the 
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counter-clockwise (^-direction. Note that at this point, indexing in 
either the radial or the (^-direction would have been equally appropri- 
ate. However, once chosen, the specified direction must be adhered to 
throughout the immediate timestep. 

(2) Two-Star Level (<j)-Direction). Because the rings asso- 
ciated with the (^-direction have neither a beginning nor an end, one 
node must be arbitrarily defined as a reference node. Using node 1 
from the one-star level and working in a counter-clockwise direction 
(Figure 8), the values of the coefficients for each node along the cho- 
sen ring were calculated and temporarily stored in suitable vectors. 
The resultant system of equations for nodes in the ^-direction was of 
the form 




Figure 8. Two-Star Level Indexing 
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(4.1) 



~b x c x 0 q- 




T,' 






a 2 b 2 C 2 0 




T 2 




' d 2 


0 b 3 c 3 




T 3 




d* 


1 

o 

i 











which is not tridiagonal. To achieve tridiagonality, the ring was cut at 
the reference node. The resulting coefficient matrix was then of the 
form 

\b x - a) Cj 0 

b 2 C 2 
0 Os b 3 

.0 0 a 4 

Node temperatures were then obtained by applying a suitable tridiago- 
nal matrix solver [Ref. 4:p. 441] and stored in the two-star intermedi- 
ate level matrix. The pointer was then incremented to the node 
adjacent to node 1 in the r-direction and the process repeated. Upon 
completion of the first layer calculations, the pointer was indexed in 
the positive z-direction and the process repeated until all two-star 
level temperatures have been calculated. 

(3) Three-Star Level (r-Direction). Construction of the 
three-star temperature matrix was accomplished in a fashion similar 
to the one- and two-star level procedures (Figure 9). Beginning at the 
reference node, the values of the node coefficients were calculated 
and temporarily stored in suitable vectors. When the coefficients for 
the last node on the selected ray had been caluclated and stored, the 
tridiagonal matrix solver was applied to the system of equations, the 



0 

0 

C 3 

C b . - c.) 



(4.2) 



38 



resultant node temperatures written to the three-star temperature 
array, the pointer indexed to the next adjacent ray in the counter- 
clockwise direction, and the process repeated until all node tempera- 
tures were calculated for the layer. The pointer was then indexed in 
the positive z-direction until all three-star level temperatures had 
been calculated. 



Z 

Q 



Y 



Figure 9. Three-Star Indexing 

(4) New Temperature Array (n+1 Level). When all inter- 
mediate matrices had been developed, new node temperatures were 
calculated in terms of the three intermediate temperature arrays and 
the array of temperatures at the previous (old) timestep in accordance 
with the discussion in paragraph 4 of Chapter III. The choice of 
indexing direction at this point was arbitrary, so the method applied 
to the star level was repeated. It must be noted that the new node 
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temperatures must be stored in a separate (new-temperature) array 
until all node temperatures have been calculated. Upon completion of 
the calculations, the new temperature array was printed or stored 
(when appropriate), copied onto the old temperature array, the time 
counter indexed, and the entire process repeated until the limit of the 
period of evaluation was reached. 

D. EXECUTION TIMEKEEPING 

The codes for both the explicit method and Brian’s Method were 
of common origin, had nearly identical main codes, and shared com- 
mon input/output and utility subroutines. The codes were organized 
such that the calculations unique to the specific method employed 
were segregated into a single, cohesive block in the form of one or 
more subroutines. A timing utility which allowed placement within the 
code of pointers to initiate and stop CPU timekeeping was employed 
to measure the execution times of each method for each example or 
application tested. Among other places, pointers were placed at the 
entry to and the exit from the respective method-specific subroutine 
blocks. Benchmark runs exclusive of I/O calls were conducted for each 
method during the validation process. Subsequent runs involving 
printed or stored output were structured such that required I/O han- 
dling instructions were common to each code to allow accurate com- 
parison of relative speeds. 
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V. VALIDATION 



A INITIAL CODE 

1. General Procedure 

To simplify code development, initial forms of the codes for 
both the explicit method and Brian’s Method were developed without 
including convection or radiation. Validation for both methods was 
accomplished by solving for the temperature distribution within the 
cylinder described in paragraph A.2 of Chapter III, comparing results 
with the analytic solution and calculating the percentage error. Upon 
achievement of satisfactory results, the initial codes were expanded to 
encompass the effects of both convection and radiation. Validation of 
the core of this revised code was accomplished by setting appropriate 
convective and radiative coefficients to zero and repeating the initial 
validation procedure with the revised code. 

2. Time Incrementation 

Initial runs were conducted using a fixed time step of 0. 1 
second. Following initial data collection, additional runs were con- 
ducted to determine operating limits for both methods. During this 
phase, Brian’s method was tested and validated using constant time 
increments in the range from 0.01 to 100 seconds and included a log- 
arithmically increasing time step arrangement summarized in Table 3. 
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TABLE 3 



VARIABLE TIME STEP EMPLOYMENT 


RANGE (SEC) 


TIME STEP (SEC) 


0 . 1 - 1. 0 


0.1 


1.0-10.0 


1.0 


10.0-100 


10.0 


100-300 


100.0 



Additional runs were conducted using a varying time step in which the 
final step size was limited to 20 seconds. 

R APPLICATION CODE 

No analytic solution was available for a cylinder with convective 
and radiative boundary conditions. Consequently, the application codes 
may not be considered fully validated. A comparison of results was 
obtained by solving the application example with both methods and 
comparing the results point-for-point against each other in the form of 
a percentage disagreement. All combinations of heat flux and boundary 
conditions considered are summarized in Table 4. 
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TABLE 4 



SUMMARY OF APPLIED BOUNDARY CONDITIONS 



Boundary 



Class 


Inner 


Outer 


Upper 


Lower 


I 


9 


R 


R 


R 


II 


9 


R9 


R 


R 


III 


9 


R9 


R9 


R 


IV 


9 


R9 


R 


R.9 


V 


9 


R9 


R9 


R.9 


VI 


9 


R9.H 


R 


R 


VII 


9 


R9 


R,H 


R 


VIII 


9 


R9 


R 


R.H 


IX 


9 


R9 


R.H 


R.H 


X 


9 


R.9.V1 


R 


R 


XI 


9 


R9 


R,V2 


R 


XII 


9 


R9 


R 


R,V3 


XIII 


9 


R9 


R,V2 


R,V3 


XIV 


9 


R9.V1 


R,V2 


R,V3 


XV 


9 


R9 


R,V2 


R,V2 


XVI 


9 


R9 


R,V3 


R,V3 



Where: 

Q = Constant, uniform surface heat flux 
R = Radiation to extensive, black surroundings 
H = Convection from surface 
V 1 * = Periodic, uniform surface heat flux 1 
V2 * = Periodic, uniform surface heat flux 2 
V3* = Periodic, uniform surface heat flux 3 



* As described in Appendix G. 
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VI. RESULTS 



A DEFINITIONS 

1. Convergence 

On validation runs, calculated values were considered to have 
converged when the net remaining change in the value comprized less 
than one percent of the analytic solution for the given time step. For 
comparison runs, values were considered to have converged when 
their difference comprised less than one percent of the value of the 
smaller of the two numbers. 

2. Steady-State 

A condition of steady-state was considered to have been 
reached when the net remaining change in a temperature comprized 
less than one percent of the final steady-state result for single tran- 
sient evaluation. 

R VALIDATION 

1. Analytic Solution 

The analytic solution converged to steady-state at 300 sec- 
onds. Consequently, all validation runs were limited to 300 seconds 
duration to economize on analysis time. 

2. Explicit Method 

The theoretical limit on the maximum time step allowed for 
the geometry selected was approximately 0.395 seconds, which did 
not lend itself to convenient comparisons with the analytic results. 
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The explicit method was validated using a numerically convenient 
time step of 0.1000 second and converged with the analytic result 
within four seconds of the start of the transient cycle (Figure 10). 
Total execution time was 40.37 seconds for a ten (radial) by four (<(>) by 
five (z) grid. The observed error at the end of the 300-second evalua- 
tion period was on the order of 0.009 percent. For processing speed 
comparison runs, the maximum theoretical timestep was employed 
and demonstrated a minimum execution time of 10.53 CPU seconds. 
Due to the odd step size, no comparison of relative accuracy was made. 

3. Brian’s Method 

Brian’s method was initially limited to a fixed timestep of 0.1 
second for the validation process to allow a relative speed comparison 
with the explicit method. While convergence times and observed 
errors were comparable (Figures 10, 11, and 12), the explicit method 
demonstrated a significant speed advantage over Brian’s Method for 
equivalent timesteps (Table 5). This situation was dramatically 
reversed when the logarithmically increasing timestep arrangement 
(Table 3) was employed and resulted in an execution speed of 2.53 
CPU seconds, which comprised a 410 percent advantage over the 
optimum result of the explicit method (Table 6). 

C FOLLOW-ON TESTING 

Upon conclusion of the validation runs and speed comparisons, a 
variety of boundary conditions with easily predictable effects were 
applied to allow observation of system response (Table 4). The follow- 
ing general observations were made. 
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10. Percentage Error of Explicit Code as Compared 
to Analytic Results and Plotted Versus Time 
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Figure 11. Percentage Error of Brian’s Code as Compared 
to Analytic Results and Plotted Versus Time 
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Figure 12. Comparison of Brian vs. Explicit Validation Results 
Expressed as a Percentage Difference 



TABLE 5 

VALIDATION RESULTS— FIXED 0.1 SEC. TIMESTEP 



10 x 4 x 5 GRID 
300 SEC TRANSIENT 



RESULT 


EXPLICIT 


BRIAN’S 


STEP SIZE 


0. 1 SEC 


0. 1 SEC 


CPU TIME 


40.37 SEC 


223.24 SEC 


ERROR AFTER 


0.1% 


15.188% 


20.308% 


1% 


1.247% 


1.668% 


2% 


0.612% 


0.814% 


10% 


0.110% 


0.183% 


100% 


0.009% 


0.009% 


CONVERGENCE 
(<1 % ERROR) 


3-4 SEC 


4-5 SEC 



SPEED ADVANTAGE: EXPLICIT (553 %) 
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TABLE 6 



VALIDATION RESULTS— LOG- VARIABLE TIMESTEP 



BRIAN’S METHOD 
VARIABLE STEP SIZE 



TIME (SEC) 

0 . 1 - 1.0 
1 . 0-10 
10-100 
100-300 



STEP SIZE (SEC) 

0.1 

1.0 

10 

20 



% TRANSIENT % ERROR % IMPROVEMENT* * 

0.1 20.308 0.00 

1.0 1.446 13.3 

10. -0.057 68.8 

100 . 0.01 - 22.2 

CPU TIME: 2.53 SEC 



SPEED ADVANTAGE : BRIAN’S (416 %) 

* Explicit method with optimum time step. 



1. Symmetric Boundary Conditions 

For both radially and axi-symmetric boundary conditions, the 
results of both methods rapidly converged with each other within the 
first three percent of the transient cycle. 

2. Asymmetric Boundary Conditions 

When asymmetric boundary conditions were applied, large 
differences were observed at surfaces nearest the surface of applica- 
tion. These differences were eliminated by apllication of a finer grid 
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size in the direction of principal flux and are thus attributed to discre- 
tization errors due to the generally coarse nature of the applied grid. 

D. APPLICATION EXAMPLE 

A satellite thermal analysis problem (Appendix G) was selected as 
offering a practical example of a system subjected to long-period 
steady-state transient behavior. Figure 13 presents a graphical 
comparison of the results of both methods for the initial 3000 seconds 
of the first transient cycle. Figure 14 exhibits the thermal response of 
selected nodes as calculated by Brian’s Method for the initial 10 tran- 
sient cycles of the application example. 
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Figure 13. Graphical Comparison of Brian’s Method and the Explicit 
for the Initial 3,000 Seconds of the First Transient Cycle of the Application ] 



SATELLITE TEMPERATURE PROFILE 

(BRIAN - 10 ORBIT) 




(o) auruvaadwai ssaoxa 
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Figure 14. Thermal Response of Selected Nodes as Calculated by Brian’s 
Method for the Initial 10 Transient Cycles of the Application Problem 



VII. CONCLUSIONS AND RECOMMENDATIONS 



A CONCLUSIONS 

1. Adaptability 

Brian’s Method is readily adaptable to cylindrical geometries. 
However, compensation must be made for the lack of tridiagonality of 
the coefficient matrix when calculating intermediate temperatures in 
the phi-direction. 

2. Relative Accuracy 

For symmetric boundary conditions, Brian’s Method delivers 
accuracy comparable to the explicit method. Convergence is negligibly 
slower than that observed for the explicit method. 

3. Relative Execution Speed 

Brian’s Method offers distinct advantages over the explicit 
method for systems demonstrating periodic steady-state temperature 
response, and for those systems with extended single transient cycles 
whose geometry makes the explicit timestep limitation overly restric- 
tive. For some geometries, Brian’s method is slower than the explicit 
method and is not recommended. 

4. Relative Benefits 

The principal benefit of Brian’s Method is its unconditional 
stability and associated freedom from timestep restrictions. This 
allows the user to apply a coarse grid and timestep to allow rapid 
coarse analysis with the ability to adjust the step size downward as the 
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situation warrants. However, the relative advantages of the unre- 
stricted timestep must be carefully considered against the substantial 
increase in algebra and code requirements. 

a RECOMMENDATIONS 

1. Adaptation of Brian’s Method FORTRAN Code to C Pro- 
gramming Language 

a Background 

The large storage and calculation requirements associ- 
ated with finite difference approximations of transient heat conduc- 
tion problems have historically relegated their analysis to the realm of 
the mainframe computer. However, recent developments in micro- 
processor technology have dramatically improved the storage capabil- 
ity and execution speeds of personal computers (PCs). Consequently, 
increasing amounts of engineering analysis and design are being con- 
ducted on these devices at significant cost reductions. In research 
being conducted under the auspices of the Electrical and Computer 
Science Department at the Naval Postgraduate School, FORTRAN 
codes involving large numbers of repetetive calculations or iterations 
have been demonstrated to experience significant improvements in 
execution speed when converted to the C programming language prior 
to execution on a personal computer [Ref. 11]. 
b. Recommendation 

It is strongly recommended that research be conducted 
to investigate the feasibility of and limitations on the conduct of tran- 
sient heat conduction analysis on personal computers employing a C 
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programming language adaptation of the Brian’s Method FORTRAN 
code. 

2. Incorporation of the Analysis of Directional and Time-Depen- 
dent Properties in Brian’s Method 

a. Background 

Composites demonstrate strong directional and 
temperature dependent thermal properties and are being rapidly 
assimilated into many engineering applications. 

b. Recommendation 

It is recommended that additional research be con- 
ducted to develop procedures for the incorporation of directional and 
time-dependent thermal properties into Brian’s Method. 
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APPENDIX A 



NODE EQUATIONS FOR EXPLICIT METHOD 
A TYPE 1 NODES— END LAYERS, INNER SURFACE 

T n+i =~^{T n +2Fo[(T'+T n ) + Cl (T n E +T H w - 2T n ) 

+ c 2 (T H 0 -T n ) 

+ c e + h c T„ + h r T sur ) 

e e 

+ C.W, +h,T.+ h, T_)J} 

where 

k A z 2 

w 

C \ A r 2 

2k z r(r+~-)&<l> 



k r (r + Af)Az* 



C 3 = 



A z 

k t 



rA z 2 

Mr + j)Ar 
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c 5 = 1 + 2Fo I c,(A + h ,) + c,(h, ' + h , ( )J 

k, =e,a(r" + T,„)[(T") 2 + T,„ r ] 

€ 

h r = £,<J (7 4- T sur )|_(7 ) +T sur J 

a TYPE 2 NODES— END LAYERS, MID-RADII 

7" +1 = + 2Fo[(T e " - 7") + c^T* + 7^ -27") 

+ c 2 (7;-7") . 

+ c ,(r;-r") 

+ c 4 W' e + h c T m + h r 7 W )J} 

where 

& A z 2 

9 

C. = 7 

2 £ 2 r 2 ^ 

k r (r- 4p)Az 2 
Ci= 2 it, rA r 2 



M^+y^z 2 

2 k z rA r 2 
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c s = \ + !Fo[c,(h,' + h,)\ 
h, = e.a(T‘ + T „)[(r - ) 2 + T J\ 

e 

C TYPE 3 NODES— END LAYERS, OUTER SURFACE 

7 " +1 =T~ 5 {T K +2Fo[{T n t -T n )+c x {T n E + T n w -2T K ) 

+ c 2 (T n I -T n ) 

+ c e + h c T„ + h r T sur ) 

€ e 

+ c.W'o+K/- + KT m )^ 

where 

U 2 2 

_ <P 

2k,(r-^f-)rA<p 2 

k r (r- ^~)Az 2 
c 2 = tr 

k : (r--^-)Ar> 
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rAz 2 



c 



4 



k,(r-^)Ar 



c 5 = 1 +2Fo 3 (h c ' + h r ) + c 4 (h c o + h r 
h r =e e o (7" + 7 w )|_(r n ) 2 +r w 2 J 
h r q = e o G (T n + T sur ) [(T n ) 2 + T sur 2 \ 

D. TYPE 4 NODES— MID-LAYERS, INNER SURFACE 

n +1 1 r n r. « n « v , n n n 

r = — {r +FoI(t n + t s -2t )+ Ci (t e + t w -2t ) 

4 

+ c 2 (r;-r‘) 

+ c 3 (<f # +A e r. + A r r w )]} 

where 

k.Az 2 

^r(r+A£)A^ 2 

2jt r (r+4r)Az 2 

C 2 = 7T 

M r + '^")Ar 2 
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2 r Az 2 



C 3 = 



k z (r+-^-)Ar 



c a =\+Fo \ c 3 (h e +h r )^ 



h r -£,<J (T +T 3 



)L o 



n ) 2 + T 5 



E. TYPE 5 NODES— MID-LAYERS, MID-RADII 

T n+1 =T n + Fo[(T n N +T n s -2T n ) + c^(T n E + T n w -2T n ) 



+ c 2 (T , - T ) 



n n I 

+ c 3 (T 0 -T )J 



where 



k. Az 2 
k z r 2 A<p 2 



M^ + y)Az 2 
k z r A r 2 



c 



3 



(r+ 4t)Az 



2 



k z A r 2 
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F. TYPE 6 NODES— MID-LAYERS, OUTER SURFACE 



7 ’ n+1 = 7 ~{T n + Fo [(r; + T] - 2 T n ) + + T n w - 27 ") 

4 

+ c 2 (T;-r") 

+ 

where 

k A z 2 

k z r(r- ^~)A(p 2 
2£,(r - ^“)Az 2 

c 2 = 77 

M r - -^p)Ar 2 

2r A z 2 

C 3 A r 

Mr- x )Ar 

C 4 = 1 +Fo |_ C 3^c 0 + /l , c )J 

hr o = e o <J (T n + T sur ) L(F n ) 2 + T SU1 . 2 \ 
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APPENDIX B 



NODE EQUATIONS FOR BRIAN’S METHOD 

A TYPE 1 NODES— END LAYERS. INNER SURFACE 

k f Az 2 

C, = -r 

2 k,r(r + -^)A <p 2 



k r (r + ^)Az 2 
k.ir+^Ar 2 



C 3 = 



Az. 

k. 



C 4 = 



rA 2 ' 



**(r + ^p)A r 



c 5 = 1 + 2Fo c 3 ( A c 



+ + c 4 (^ C/ + 



h r = e e o{T n + T s j[(T n ) 2 + T S J] 
h r =£,G(T n +T s j[(T n ) 2 + Tj] 
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1. Step 1 (z-Direction) 



AT\+BT m +CT' n =D 



where 



A = 




Bottom layer nodes 
Top layer nodes 



B 



1 + Fol \ + c 2 (h Ct 




C = 




Top layer nodes 
Bottom layer nodes 



D = T K + Fo[c l (T K E + T H w -2T' , )+c 2 (T n 0 -T 

+ C -}W' e + h c T„ + h r T sur ) 
e t 

+ c,w',+h,T. + h r r_)J 

2. Step 2 (^-Direction) 

at" +bt" + ct” =D 



where 



A = - c, Fo 
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B = \+Fo ]2 c, + c 2 (h C ' + h r )+c A (h C/ +h r/ )J 



C ■= — c x Fo 



D = T + Fo La". - T m ) + c A (T n 0 - T n ) 



+ c 3 (</ / , + h c T „ + h r T lur ) 

+ c,(<f,+h,T. + h, r„)J 

3. Step 3 (r-Direction) 

AT +BT“' + CT“‘ =D 

where 



A= 0 



B = l + Fo 



c +c (h c +h r )+c A (h c +h r )\ 

_ e e I 1 J 



C = - c 2 Fo 



d = t n + Fo [(r; - r* ) + c, (r;‘ + r*; - 2r“) 
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+ cj<r,+h'T. + h, T„)j 



4. Step 4 (^-Direction) 

r" +1 =-^-{r"+2Fo[(r* < + r*) + c 1 (7'y + r‘;-2r") 

+ c 2 (r;--r**) 

+ c 3 (</' , + /z c r» + /i r F w ) 

< < 

+ C.WV*. T.-h, r„)J} 

a TYPE 2 NODES— END LAYERS. MID-RADII 

k A z 2 

c i = 7 

2k z r 2 (p 2 

Ur 

2 2k z rAr 2 

kXr + ^)Az 2 

C, = ; 

3 2k z rAr 2 

A z 



c 5 =l+ 2 Fo[c 4 (/ l£< + /z r< )J 
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h, = £ ,o <r +7-„)L<rV + r, 2 ,J 

€ 

1. Step 1 (z-Direction) 

AT’ S +BT' +CT' n 



where 



A = 




Bottom layer nodes 
Top layer nodes 



B= 1 + 



Fo[1+ c 4 (/? Ce + 




C = 




Top layer nodes 
Bottom layer nodes 



D = T n + Fo[c : (T n E +T n w - 2T n ) + c 2 (T n , - T n ) 



+ c,(T 0 -T n ) 



+ c t + h e T„ + h r T w )J 



2. Step 2 ((^-Direction) 



at'* +5r’* + cr“ =d 



where 



A = - c, Fo 
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B — 1 + Fo |^2 c , + c 4 ( /z c ^ + /i r ^ ) J 



C = - c x Fo 



D = T' + Fo[(T;-T')+ c,(T" -T') + c ,( T" 0 - T ■) 
+ c < (/. + A I< T. + A,J„)J 



3. Step 3 (r-Direction) 



at, +sr +cr 0 =£> 



where 



>4 = - c 2 Fo 



B= \+Fo^c 2 + c 3 + c A (hc e + h r)\ 



C = - c 3 Fo 



D = T " + Fo [a; - T ‘ ) + c , (r ** + T ” - IT ** ) 



+ c 4 (^' t + /i £ + h 



r 
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4. 



Step 4 (New Temperatures) 



T 



n+1 



= {t n + 2Fo [(t: - t ) + c .(r;* + r; - it ** 

+ c 2 (t‘* - r‘” ) + c 3 (r“* - r’** ) 

+ c A W' t + h c T m + h r r w )J} 



C TYPE 3 NODES— END LAYERS. OUTER SURFACE 

& a A z 2 

2A: 2 (r-^)rA^ 2 



M r - -4p)Az 2 
Mr--^)A r 2 



c 



3 



Az 



C 4 = 



rAz 2 



t A r. 
,('-T)Ar 



c 5 = 1 + 2F<?| c 3 (/j + A r )+c 4 (/» c +h r )| 

L « < o o J 



*, = e,c(r*+ r,„)[(TV + r,„ J ] 
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h, 0 = £ 0 o (T n +T sur )[(T n ) 2 + T sur 2 \ 

1. Step 1 (z-Direction) 

AT] +BT' + CT] =D 



where 



f 0, Bottom layer nodes 

A. ~~ i 

Fo, Top layer nodes 

B = l + Fo\l + c 3 (h c +h r )+e A (h c +h r )l 

- Fo, Top layer nodes 
0, Bottom layer nodes 




D = T~ + Fo[c^T' c +T'„ -2T") + c 2 (T; -T') 
+ c,(<f. + h, T.+ h 7 T ) 

e e 

+ c A (cf Q - h c T„- h r T sur )J 

2. Step 2 ((^-Direction) 

AT]‘ +5r“ + C7’’,‘ =D 

where 

A - - Cj Fo 
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B = 1 + Fo |2c, + c 3 (/i c< + h r ) + c 4 (h + h r J J 



C = - c x Fo 



D = T" + Fo[(Tl-T m ) + c 2 (Tj -T”) 



+ Cj(<f, + h c T„ + h r T„) 

€ e 

+ o + h + h r T sur )J 



3. Step 3 (r-Direction) 



at]" +bt'" + ct"‘ =d 



where 



A~ - c , Fo 



B = 



1 + Fo |^c 2 + c 2 (h c * + h r ) + c 4 (h c o + h T J j 



C = 0 



D = T n + Fo [(T‘ -T‘) + c x (T" + T"-2T") 



+ c i(<f' t + h c T„ + h r T sur ) 
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+ c M' O +h c 0 T - + h r T sur ) J 



4. Step 4 (New Temperatures) 

r" +1 =~t s {t k +2Fo[(t]~ r‘) +c 1 (T‘* + r;*-2r'*) 



+ c 2 (r;**-T“*) 

+ c 3 (</' « + ^c T„ + h r T tur ) 

t t 

+ e + h c T „+ h r J sur )Jj" 

D. TYPE 4 NODES— MID-LAYERS, INNER SURFACE 

k Az 2 
_ * 

k t r(r+4£) A<P 2 



2 k r (r + -4p)Az 2 
k t {r + -^p)A r 2 



c 



3 



2 rA z 2 

Mr + — ^A r 



c 6 = 1 +Fo [c 3 (h c ^ + h r )] 



h 



= £,C ( 7 + T a 



■)Lo 



v + t 3 
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1. Step 1 (z-Direction) 



AT] +BT' +CT' n =D 



where 



A = 



Fo 

2 



B = l + f %-[2 + c 3 (h Ci + h r )] 



C 



Fo 

2 



D = T n + f f L [c 1 (T n E+ T n w -2T n ) + c 2 (T n 0 -T n ) 



+ c i i<f,+h e T m + h r T m „) J 



2. Step 2 (^-direction) 



AT‘‘ + 57'** + CT~ =D 



where 



A 



c, Fo 

~T~ 



fi = l + “[2c 1 + c 3 (* e/ +A r/ )] 
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c = - 



c, Fo 
~2 



D=T' +^-[a' rl + T' s -2T')+c 2 (T^-T") 
+ €,(.</' ,+h,T. + h,T„) J 



3. Step 3 (r-Direction) 



AT]" +BT'" + CT *“ = £> 



where 



A=0 



Fo r 

2 L “ 2 ' " 3 



B = 1 + -T- c, + c 



(A e/ + /» r/ )] 



C = — 



c„ Fo 



O = r + -^ [<r ; + r; - 2T • ) + c , (T " + T " - 2T ") 



+ c ,(<r ,+h.T. + h r r_)J 



4. Step 4 (New Temperatures) 

r " +1 = 77 { r " + + t' s - it * ) +c,(r;’ + r;* - 2r *’) 
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* • • • • * , 



+ c 2 (T 0 -T ) 



+ C,W , + KT- + h,T„) J} 



E. TYPE 4 NODES— MID-LAYERS, MID-RADII 

A z 2 

C ' k t r 2 y > 2 

k r (r + 4p)Az 2 
° 2 k t rAr 2 



C 3 = 



(r + -^)Az 2 



k t Ar‘ 



1. Step 1 (z-Direction) 

AT] +BT‘ +CT' n =D 

where 




B = 1 + Fo 



C = - 



Fo 

2 
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P f) r n n n n n nn~l 

D = T + -Y [ Cl (T E + T w - 27 ) + c 2 (T i - 7 ) + c 3 (7 0 -7 )J 



2 . 



where 



3. 



where 



Step 2 (^-Direction) 

at e + 57“ + C7~ =D 



A = 



Cj fb 

~T~ 



5=1+ c^Fo 



C = - 



c, fb 

"T - 



o = r + -^[<r; + T;-27-') + c 2 (r;-r") + c J (r;-T")] 

Step 3 (r-Direction) 



47, +57 + CT „ =D 



A = 



c 2 7o 
2 



5 = 1 + ^ (c 2 +c 3 ) 
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c = - 



c 3 Fo 
~ 2 ~ 



D = T’ , + f j L [(T; + T' s -2T') + c l (T' E ' + T';-2T‘‘)] 



4. Step 4 (New Temperatures) 

7’" +I =T " + Fo (V; + T' S -2T') + Cl (T‘‘ + T” - 2 T **) + c 2 (7’ / ‘“ - T *“) 



+ 




F. TYPE 4 NODES— MID-LAYERS, OUTER SURFACE 

k $ Az 2 

k, r(r- -^)A<p 2 
2k r (r- ~~)Az 2 

c 2 a7 

^(r--j)Ar 2 

2 rA z 2 

C 3 — a r 

k t { r — — )A r 

c 4 = 1 +Fo|^c 3 (/ 2 Co + /i ro )J 

hr Q = e 0 o (7' n +7’ Xkr )l_(r n ) 2 + r w 2 J 
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1. Step 1 (z-Direction) 



where 



AT’ S +BT +ct' n =d 



A = - 



Fo 

2 



5=1 + 





+ c 2 ((/' o + ft ‘ 0 T “ + hr 0 T ^ 

2. Step 2 (^-direction) 



where 



AT” +BT “ + CT” =D 



A 



c, Fo 

~T~ 



Fo 



B = 1 + — + 2c, + c 3 (h, 



+ h r ) 

O O 



c = - 



Cj Fo 

~T~ 
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D = r" + -y [<r; +r;- it ') + c 2 < t "- t ") 

+ c,(<f 0 + h t T.+ h r T ) I 

o o J 



3. Step 3 (r-Direction) 



AT, +BT 



+ CT 0 =D 



where 



r Fn 




5 — 1+ c + c Ah c +h r ) 
l 1 o o 



C = 0 



D = T 



+ 



-T-[(T' n +T' s - 2T') + C] (T n E +T n w -2T n ) 
+ c 3 i<f 0 +h e T m + h r T„) J 



4. Step 4 (New Temperatures) 



T"' =j-{T- + F„[(,T;+Tl-2T) + c l (T-; + T-;-2T‘ 



+ c 2 ( 7 / -T ) 
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APPENDIX C 



VALIDATION MODEL 
A PROBLEM SPECIFICATION 

Given a right circular cylinder with adiabatic ends and specified 
initial temperature distribution and constant inner-surface tempera- 
ture of zero degrees Celsius. At time t = 0, a uniform constant heat flux 
is applied to the outer surface of the cylinder. Assume constant prop- 
erties and no losses due to convection or radiation. 

B. ANALYTIC SOLUTION 
Fourier Equation 



vY = -1- — 

« a 



Boundary Conditions: 



r = a 



T = 0 



r D = b 



dt k 



Initial Condition : 



t = 0 



7 = 0 



Let: T= r„ t) + <*>( r c ) 
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Auxiliary Problem 



, , i dWr-.t) 

V" f'C r 0 ,t) + V‘cP( r 0 ) = -g 



r = a 



•Kr 0 ,t) + V 2 <*>(r 0 ) = -^^ 



r n —b 



& (r..O + #-=-T 

dr K 0 ' dr 0 k 



t = 0 



•F ( r o ,0)+® (r,) = 0 



Separating the Auxiliary Problem 
Equation 1: 



d® , X n 
r ' 



r -= a 



0 (a) = 0 



r,= b 



Let: 



y = 



d® 
dr o 



Then: 



d® _ q" 

dr B k 



dy 

dr o 



y 
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dy 



y ~ 



dr 

w o 



In ( y) =-ln(r c ) + c, 







And: 



dO c i 

dr 0 ~ C 



d& = Cj 



dr 

o 



r 



o 



Yielding: 



d> =c j ln(r e ) + c 2 



Evaluating Constants 



r = a -» <2> =0 



0 = Cj ln(a)+ c 2 



C 2 = ~ 



* 



ln(fl) 



dd> _ Q" 
dr a k 



b d' r 

Thus: 0(r .)=-%- In (-^) 



c i = 



Jfc 
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Transforming the Variable 

Let: r=-x 

b 

And: K= ^ 



btf r 

Then: 



Converting Boundary Conditions 



r o = a => 




K 



r »- b => r = 1 
The Transformed Problem Becomes 

bcf r 

0(r)= ~r 



r — K <P = 0 

d<K r) bq" 

k 

1 dH' 1 dV 

dr, 1 + r ° dr B « dt 

r c = a na)= 0 



dr 

Equation 2: 

<?¥ 
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r,= b 



dV_ 

fr 0 



= 0 



t = 0 (r o ,0) =- O (r e ) 



Let: <F ( r o ,t)=R(r 0 )r (: ) 

Then: ^ F + y^r=^-^r 

Which reduces to 



M / 

R r o R 



_l /i 
a r 



= -P * 



Separating the two equalities, we obtain 
Equation 2. a 



L iL 
a r ~ 



P 



. 2 



In (r )=- a P ,2 t 



r«) = e- a 



p 



*2 



l 



Equation 2.b 

r o= a R = 0 




From Table 3-2, Case B of Ozisik 
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R(P \r 0 )=JAP * r 0 ) Y' 0 {p ' b) - J AP * b) ■ YAP ’ r c ) 



the inverse of the norm is given by 

j ,q) 

N (P *) ~ Jl(P‘ a)-j 2 0 (p -b) 
and the eigenvalues are the roots of 

JAP * a)-Y 0 (p * b)-J 0 {p'b)-YAP *a) = 0 
From the recursion relations we know 

J' c (x)=- J{x) 

and similarly 

FlU)=- F,(jc) 

Substituting, multiplying by -1, and rearranging terms yields 
y, (p * b) J a (P * a)-J x {P -b) Y 0 (P -a) = 0 

Transforming Variables: 

Letting: P = P b 
and recalling that: k = 

we obtain by substitution the common tabular form: 

Y ] (P) JAk p)-JAP )-y,(k p ) = 0 
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Performing similar variable transformations on the transient 
solution yields: 



RW,r) = J l (p)-Y{rP)-J 0 (rp)-Y l (p ) 



with the norm 



* JUv/n 



Total Solution: 

b q r t-> — oc B 2 t i 

T(r,, ) = h (i ) S (exp (— ^-) * . i'.P ) 

• Ap *,(>%/? ) In (f )] rfp } 
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APPENDIX D 



EXPLICIT METHOD STABILITY REQUIREMENTS 

The following equations represent the limiting time step size for 
their respective node types. 

A TYPE 1 NODES— END LAYER, INNER SURFACE 

k A 2 2 k r (r +Ar) A z 2 

A r, <1 + — - + j 

k ,r(r + -^)A<p 2 k J (r + -^-)Ar 2 



B TYPE 2 NODES— END LAYERS, MID-RADII 



k A z 2 k r (r - Af~) A z 2 

Af < 1 + — T ± = 

k, r 2 A(p 2k,r A r 



C TYPE 3 NODES— END LAYERS, OUTER SURFACE 



k A z 2 k r (r + Af-) A z 2 

A r 3 < 1 + *— + f 

k z r(r + -^)A(t > 2 k t (r+^)Ar 2 



D. TYPE 4 NODES— MID-LAYERS, INNER SURFACE 



2 k A z 2 2 k (r + ) 

At, < 2 + + — 

k'(r--^)A<t> 2 Ar 
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E. TYPE 5 NODES— MID-LAYERS, MID-RADII 



2k A z 2 k r (r - A~) k r (r+ A p) 

<2 + - =- + 4 — + 4 — 

k t r 2 A(p k t rAr 2 k t rAr 2 



F. TYPE 6 NODES— MID-LAYERS, OUTER-SURFACE 






2k f Az 2 



2k r (r- A £-)Az 2 



k t r(r-^) Atp 2 t,(r- j)Ar : 
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APPENDIX E 



c 

c 

c 

c 

c 

c 

c 

c 



c 



c*** 



c 



c 



c 



c 



EXPLICIT METHOD CODE 



ECR FORTRAN A - (REV. 11/23/88) 

SATELLITE MODEL - VALIDATED 11/26/88 
********************************************************** 

* ECR - PROGRAM TO COMPUTE NODE TEMPERATURES FOR A RIGHT * 

* CIRCULAR CYLINDER SUBJECTED TO CONVECTIVE AND/OR * 

* RADIATIVE BOUNDARY CONDITIONS USING THE EXPLICIT * 

* METHOD OF FINITE DIFFERENCE ANALYSIS. * 

********************************************************** 

IMPLICIT DOUBLE PRECISION (A-H,0-Z) 

DIMENSION TOLD (15, 15, 15) ,TNEW(15, 15, 15) 

DIMENSION DMAT (15, 15, 15) 



COMMON /ARRSIZ/ III,JJJ,KKK 
COMMON /CONDUC/ DKR, DKP, DKZ 
COMMON / CONVEC/ HSUBN, HSUBS , HSUBI , HSUBO 

COMMON /GEOMET/ DPHI , DR, DZ , RAD, RMAX, RMIN, THICK, ZMAX, LPHI 
COMMON /HEAT/ ALFA, FO, QPSUBN, QPSUBS, QPSUBI, QPSUBO 
COMMON /MATSIZ/ NODES , NPHI , NLAYER 
COMMON /PARAMS/ CSUBP,RHO 

COMMON /RADN/ EPSUBN, EPSUBS , EPSUBI , EPSUBO, S IGMA 

T1,T2,T3,T4,T5,T6, T7,T8,T9,T10, TINF, TSUR 
DELTIM, DTMAX, FINTIM, LIMIT 

START TIMER 
CALL STIME 

OPEN (UNIT=1 1 , FILE= ' ECRLNG ' ) 

OPEN (UNIT=14 , FILE= ' ECRDAT ' ) 

OPEN (UNIT=15 , FILE= ' ECRCPU ' ) 

******************** USER DEFINED VARIABLES ************************ 



COMMON /TEMPS/ 
COMMON /TIMES/ 



DELTIM = 0 .30D+00 
FINTIM = 3000 . 0D+00 
IPRINT = 100 

ACTUAL MATRIX DIMENSIONS 
NODES =10 
NPHI = 4 
NLAYER = 5 



NODE SPECIFICATION FOR FILE ESHORT 
MRAD = 10 
MPHI = 2 
MLAYER = 3 

DECLARED ARRAY DIMENSIONS 
III = 15 
JJJ = 15 
KKK = 15 

ABSORPTIVITY COEFFICIENTS 
ALFAN = 0 . 4D+00 
ALFAS = 0.4D+00 
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ALFAO = 0.4D+00 
C ALBEDO COEFFICIENT 

ALBEDO = 0.30D+00 
C THERMAL CONDUCTIVITIES 

DKR = 177.D+00 
DKP = 177.D+00 
DKZ = 177.D+00 
C EMI SSIVITIES 

EPSUBN = 0.9D+00 
EPSUBS = 0.9D+00 
EPSUBI = 0.D+00 
EPSUBO = 0.9D+00 
C CONVECTION COEFFICIENTS 

HSUBN = 0.D+00 
HSUBS = 0.D+00 
HSUBI = 0.D+00 
HSUBO = 0.D+00 
C HEAT FLUXES 

SOLAR = 1353. D+00 
EIR = 238.0D+00 
GEN = 1000. 

QPSUBI = 0 .D+00 
QPSUBO = 0 .D+00 
QPSUBN = 0 .D+00 
QPSUBS = 0 .D+00 
C GEOMETRY 

RMIN = 0 . 15244D+00 
RMAX = 0 . 228 66D+00 
ZMAX = 1 . D+00 

C VIEW FACTORS WITH THE EARTH 

VERTHN = 0 .D+00 
VERTHS = 1 .D+00 
VERTHI = 0 . D+00 
VERTHO = 0 .D+00 
C VIEW FACTORS WITH SPACE 

VSPCN = 1 .D+00 
VSPCS = 0 .D+00 
VSPCI = 0.D+00 
VSPCO = 1 .D+00 
C TEMPERATURES 

TABS = 273. D+00 
TINIT = 293. D+00 
TINF = 0 .D+00 
TSUR = 4. D+00 
C MATERIAL PROPERTIES 

CSUBP = 875.0D+00 
RHO = 2270 . 0D + 00 
C INITIALIZE OTHER VARIABLES 



I COUNT 


= 


0 


ITIME 


= 


0 


ORBTIM 


= 


0 .D+00 


PI 




4 . 0D+00*DTAN (1 .D+00) 


RAD 


= 


RMIN 


SIGMA 


= 


5.67D-08 


TIME 




0 .D+00 
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T1 = O.D+OO 

C INITIALIZE MATRICES 

CALL BMINIT (TOLD, TINIT) 

CALL BMINIT (TNEW, TINIT) 

C SURFACE AREA CALCULATIONS 

SAREAN = PI* (RMAX**2 - RMIN**2) 

SAREAS = PI* (RMAX**2 - RMIN**2) 

SAREAI = PI* (2*RMIN) *ZMAX 
SAREAO = PI* (2*RMAX) *ZMAX 
C INITIAL CALCULATIONS 

ALFA = DKZ/ (RHO*CSUBP) 

DPHI = 2 . *PI/NPHI 

DR = (RMAX-RMIN) / (NODES-1) 

DZ = ZMAX/ (NLAYER-1) 

LPHI = 360/NPHI 
QPSUBI = GEN/ SAREAI 
C ORBITAL PARAMETERS 

PERIOD = 5440. D+00 
ECLIPS = 2181 . 9D+00 
OMEGA = 1.155D-03 
TAU = 0 . 01*PERIOD/4 
TIME1 = 0.D+00 
TIME2 = PERIOD/4 .D+00 
TIME3 = 107 . 8D+00/360 * PERIOD 
TIME4 = (270 - 17.8) /360 * PERIOD - 4*TAU 
TIME5 = (270 - 17 . 8) /360 * PERIOD 
TIME6 = 270/360 * PERIOD 
C ESTABLISH STABILITY REQUIREMENT 

CALL STABIL (DTMAX) 

FO = ALFA*DELTIM/ (DZ**2) 

LIMIT = FINTIM/DELTIM + 500 
C WRITE MATRIX DIMENSIONS TO FILE 14 (ECRDAT) 

WRITE (14,*) NODES, NPHI,NLAYER 

*** PRINT OPTION 3 - SAVES DESIGNATED NODE TEMPS AT SPECIFIED 

PRINT INTERVAL. USE FOR LONG PERIOD EVALS . 

WRITE (11,*) ' ECRSAT ' 

WRITE (11,*) 

WRITE (11,888) 

888 FORMAT (T6, ' TIME ' , T18 , ' DELTIM ’ , T30, ' LAYER1 ' , T45 , ' LAYER3 ' , T 60 , ' LAYE 
1R' ) 

WRITE (11,*) 

C 

CMMMMMMMMMMMMMMMMMMM BEGIN TEMPERATURE COMPUTATIONS MMMMMMMMMMMMMMMMMMMM 
C GET TIME ENTERING LOOP (MAINFRAME VERSIONS OF CODE ONLY) 

WRITE (11,*) 'ENTERING LOOP TIME' 

CALL GTIME ( 1, 11) 

DO 100 1=1, LIMIT 

ICOUNT=ICOUNT+l 

TIME=TIME+DELTIM 

ORBTIM = ORBTIM + DELTIM 

IF (ORBTIM. GE. PERIOD) ORBTIM = ORBTIM - PERIOD 
C SURFACE FLUX CALCULATIONS (TIME DEPENDENT) 

IF ( (ORBTIM. GE.0) .AND. (ORBTIM. LE.TIME2) ) THEN 
QPSUBN = ALFAN*SOLAR*DCOS (OMEGA*ORBTIM) 
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